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RESEARCH  OBJECTIVE 


The  research  performed,  and  reported  on  herein,  is  a 
continuation  of  research  performed  at  Rome  Air  Development 
Center,  Griffiss  AFB,  during  the  Summer  of  1981  under  the 
USAF-SCEEE  Summer  Faculty  Research  Program  |1|.  This 
research  is  therefore  directed  at  the  enhancement  of 
resolution  in  scenes  with  limited  support.  The-  importance 
to  the  Air  Force  lies  in  the  application  to  space  based 
infrared  sensors.  In  fhi..  context  we  have  a  field  ol  view 
limited  by  a  sunshade,  operated  on  by  a  Fourier  transforming 
lens,  and  subsequently  sampled  over  a  finite  support  by  a 
detector  array. 

The  idea  is  to  reconstruct  the  limited  field  of  view 
scene,  most  compatible  with  the  given  information.  This 
information  consists  of  frequency  domain  samples,  ultimately 
arrived  at  by  measurements  and  therefore  corrupted  by  noise, 
in  addition  to  a  priori  information.  The  a  priori 
information  yields  constraints  on  solutions  to  the  problem, 
by  such  requirements  as  the  nonnegativity  of  the  scene 
function,  the  known  physical  extent  ol  the  scene  function, 
and  probabilistic  characterisations  of  the  corrupting  noise. 

The  objective  of  the  proposed  research  then,  is  to 
evaluate  the  performance  of  iterative  deconvolution 


algorithms,  in  the  context  of,  and  with  the  constraints  for, 


the  space  based  infrared  sensor  application.  The 
sensitivity  of  the  solutions  with  respect  to  the  various 
constraints  will  indicate  the  robustness  of  the  algorithms 
against  assumptions  made. 

Until  now  iterative  reconstruction  algorithms  did  not 
incorporate  effective  ways  of  dealing  with  measurement 
noise.  Consequently,  the  noisy  measurements  obtained  in  a 
practical  application,  could  be  incompatible  with  the  a 
priori  constraints.  This  then  leads  to  searching  for  a 
nonexisting  solution  to  the  formulated  problem.  In  this 
stage  of  the  research  we  concentrate,  therefore,  on 
developing  modifications  of  the  original  algorithm  that 
explicitly  allow  the  use  of  knowledge  about  the  noise 
corrupting  the  measurements.  Analysis  of  the  resulting 
"soft  constraint"  algorithm  shows  that  the  essential 
convergence  properties  are  preserved.  initial  simulation 
results  show  the  robustness  of  the  new  algorithm,  when 
compared  to  the  traditional  "hard  constraint"  version. 
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RESEARCH  STATUS 


1 .  INTRODUCTION 

In  our  problem  an  optical  sensor  has  a  sunshade  that 
limits  the  field  of  view.  The  limited  scene  is  subsequently 
Fourier  transformed  and  a  number  of  samples  is  taken  in  the 
Fourier  domain.  These  frequency  domain  measurements:  c ii>- 
naturally  subject  to  measurement  noise.  The  idea  now  . ..  to 
recover  as  much  as  possible  of  the  original  scene. 

Philosophically  at  least  this  is  a  re  a.  onui.  it-  ;  den  due  ’ 
the  analogy  with  the  problem  of  extrapolating  ba.-.u. .  m.  t 
signals.  Several  iterative  projection  algorithm-  };«•.••..  been 
proposed  |2],  and  proofs  of  convergence  have  been  giver,  for 
the  continuous  problem  | 3 , 4 | .  The  successful  extrapolation 
by  a  factor  of  20,  in  (5  1,  led  the  author  to  extend  and 
apply  that  computationally  efficient  one-shot  approach  to 
the  2-D  problem  outlined  above  J  1  )  .  The  results  of  that 
approach  were  not  encouraging  due  to  an  extreme  sensitivity 
to  noise,  as  well  as  the  difficulty  of  implementing  all  a 
priori  information.  An  iterative-  projection  algorithm  (A! 
facilitates  the  incorporation  of  a  priori  information  at  the 
expense  of  an  increased  computational  burden.  We  develop 
a  modification  of  some  of  the  constraints,  to  accomodate 
the  effects  of  measurement  noise. 
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I 


2.  PROBLEM  FORMULATION 

As  extensions  to  L-D  are  rather  elementary,  1-D  formula¬ 
tions  are  used  for  transparency  of  the  derivations. 

Let  f(x)  denote  the  continuous  space-limited  scene,  and 

F  (u)  its  Fourier  transform.  The  sequence  {F  }  is  obtained 
c  ‘  n 


from  F^lw)  by  sampling  with  interval 


A  periodic  function 


f  (t),  period  2  n,  can  be  associated  with  {F  }  as  follows 
s  n 


f  (t)  =  Z  F  e 
s  n 

n —  °° 


•  jtn 


(1) 


For  the  Fourier  cn 


W  *  *  1  *  .  1  .  t  ' 


F  =  — 
n  2  n 


f  ( t )  e  ^  "ndt 


(2) 


-a 


As  {F  }  was  obtained  by  sampling,  the  following  relationship 


is  valid 


F  =  F  ( n2 ) 
n  c 


f  (x)  e^m‘Xdx 


<»  (  2  r  +  1  )  :  /\i 

v  /  f(x)e^nAXdx 

r=—  (  2  r  -  1 )  '/SI 


Substituting  x 


for  x  y  1  e  .1.  d  s 


F  -  r 
n 


J 


t  ,.-,221,  -indx  -  ] ndd  r 
I  f(x  +  -j7*>e  e  dx 


r —  —  -it/2 

Recognizing  the  last  exponential  to  have  a  power  which  is  an 
integer  multiple  of  2  rr,  and  substituting  t  for  xS2,  results 


in : 


f ( ( t+2rr )/2)  e~]nt  dt/fi 


(3) 


r--oc 


•  i 


Comparing  (2)  and  (3)  yields  the  relationship  between  the 
scene-limited  function  f(.)  and  the  periodic  function  f r  (  .  ) . 


f  (t)  =  —  ::  f  ( ( t+2  r)A) 

r=-~ 

The  corresponding  2-D  result  follows. 


(4) 


fs(tl’t2) 


r2= 


,  ^  ( ( t  L+2  r L , 


(t2+2"r  V:.2) 


Assume  that  f(x)  is  scene- 1 imi ted  to  ix.,x0|,  an  interval 
length  x2~x^.  From  (4)  the  original  function  f(x)  can  oe 
recovered  from  the  periodic  function  f s ( t )  associated  with 
the  frequency  domain  samples,  if  those  samples  are  spaced 
closely  enough,  that  is 


The  important  result  is,  that  the  space- 1 imi ted  scene  can  be 
recovered  from  all  frequency  domain  samples,  under  the 
constraint  of  (5).  In  the  sequel  (5)  will  be  assumed  valid 
by  some  margin,  so  called  oversampling,  so  that  fo(t)  is 
truly  scene- 1 imi ted  to  a  subset  of  [-  ,  ■]. 

It  is  a  real  world  problem  that  all  frequency  domain  sam¬ 
ples  are  necessary.  For  the  space  based  infrared  sensor 
this  requirement  translates  into  an  infinitely  large  detec¬ 
tor  array,  which  is  clearly  out  of  the  question.  In  prac¬ 
tice,  therefore,  we  hope  to  get  access  to  the  dominant  portion 


of  iF  ),  by  means  of  an  extrapolation  process,  so  that 


-n 


n=-o 


N 


-  n 


A  uniform  sampling  of  f  (t)  on  the  unit  circle  in  the 
plex  i-plane,  yields  an  N-sequence,  { f  }  say,  for  which 


{(>) 


fk  =  fs(T)/i  =  w"k 

N 


=  y  v  w 


xn 


n=-» 


n. , 
n  N 


where 


w  £  e-j2/N 

N  e 

Now  a  1-to-l  relation  exists  between  the  N-sequence  1 f 

'  i 

the  N-sequence  { F^, 'l . 

,  N- 1 

c  -  A.  V  f  \-J  ~  n 
n  N  ^._-q  x  N 

A  substitution  of  (3)  into  (  10)  yield:, 

N-l  « 


1_ 

N 


k  =  0  r.~- 

r  .  N-l 


,,  ,  .k:..,  ,-xn 
1  W  w 
rn  N  N 


m 


k~0 


-k(n-m) 

N 


(«) 


(  v  l 


ana 


(10) 


m=-» 


m  m-n-rN 


V  integer  r 


r — -oo 


’  n  +  rN 


an 


Consequently  (7)  will  be  an  increasingly  better  approxima¬ 
tion  as  N  becomes  large  enough  for  the  r-0  term  in  (11)  to 
dominate.  The  accuracy  of  the  detail  in  the  reconstructed 
scene  depends  on  N  and  the  degree  to  which  (7)  is  satisfied. 
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A  graphical  repr  i  sent  a*  : cn  -  :  problem 

Figure  1 . 


3.  SOLUTION  AFPROACH 
A.  Alternating  F  rejection  Algorithm 

For  the  derivation  and  anaiys;.  that  1  •  ■ :  :  <n  v.r :  ,r  v. 

(6  1)  we  assume  the  discrete  frequency  sequence  <  F  }  a.nci  in- 
asscciatea  periodic  function  f  (t)  over  \-r,-  I  : crnstitute 
our  problem.  Only  a  cm. ill  numbei  c:  trie  i  r r.  .y  .:oir«ir: 


samples  is 

available  due  tc  , 

;  lr  - 

!  i  - ' )  1  ic  n*  <  i  i  r.  •  i . 

’  ■"  -  c 

truncation 

ODorator  7 . 

V  -  T  F 


(  ) 


The  aim  new,  is  to  use  all  available  a  priori  iniormation  to 
extrapolate,  or  rather  estimate,  the  i  requeue  y  domain  cam¬ 


ples.  If  such,  ent  i:aat  ion  l .  r.uvoirc.fu  l  ,  a  1  •  t  *  *  •  i  •  •stnrm’  •• 
of  the  original  is  achieved  than  would  be  possible  from  the 
measurements  without  extrapolating  estimation. 

In  addition  to  (If)  a  solutioi.  F  must  satisfy  a  number  of 
constraints  ir.  the  scene  domain,  indicated  by  n  compound  op¬ 
erator  C  ■  Also  a  constraint  in  the  freouenev  domain 
s 

ists,  and  is  denoted  by  C„.  Consequently,  we  have 


F  =  FC  F~ 1 C  F 
s  F 


(13) 


or  equivalently 


F  =  c  f c  r-F 

F  s 


(  14) 
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For  rotational  convenient---; ,  -t-t'..  define  a  c  n.  tin.nc 

tor  f  encompassing  (13)  and  (14),  so  that 

F  =  OF  (13) 

From  (L2)  and  (15)  the  following  equality  is  valid,  lor  any 

F  =  CF  +  u  (Y-TCF) 

=  OF  (16) 

=  pY  +  (  i-v  T)  CF  (It) 

The  problem  is  to  find  a  solution  F  to  the  abo-e,  or  in  oth¬ 
er  words,  to  find  a  so-called  fixed  point  of  the  operator  0 . 

A  common  mode  of  attack  is  the  successive  approximation 
approach,  implemented  by 

Fk+i  =  uY  ♦  (I-,T)0Fk  (16) 


According  to  the  contraction  mapping  theorem,  (18)  yields  a 
unique  solution  if  t1  is  a  contraction  operator,  i.e.,  if 
0  c.  <  a.  i  exists,  such  that 


'r0F2;  i<  tl  !FrF2i  iVFlfF2 


(  19) 


If  a  can  take  on  the  unit  value,  (14)  d«.-i  men  a  nonexpans i ve 
operator.  A  nonexpansive  operator  there! ore  decreases  the 
distance  between  signals.  The  equivalent  requirement  for 


convergence  of  (13),  is  that  (  I-v.  7  )d  a  (ontracrtion  opera- 


tor.  For  the  latte*  ;  *.  : s  .  ..  i.-n:  t  hut  1  - 1  ana  C  \y- 

both  nonexpans  1  ve ,  and  one  ot  t  i  va  contraction  opera¬ 
tor.  As  the  second  term  or.  the  riyhthand  side  of  (16) 

constitutes  a  correction  term  over  the  measurement  domain, 
one  often  encounters  the  choice  y=l,  expressing  no  need  for 
correcting  yY  over  that  domain.  We'll  see  shortly  that  in 
the  present  approach  such  a  choice  is  not  applicable. 

B.  Soft  Constraints 

One  of  the  most  crucial  issues  in  having  a  chance  at  es¬ 
tablishing  some  kind  of  convergence  for  the  iterative  pro,. 
in  (18),  is  the  compatibility  of  the  r,<-t  of  constraints  that 
applies  to  the  solution.  As  physical  measurements  are  inva¬ 
riably  subject  to  noise  effects,  it  seems  ill-advised  to  im¬ 
pose  these  noisy  measurements  as  a  hard  constraint,  i.e.  as 
if  these  constituted  absolute  kn>  w  i  ■  ■  ig«  • .  As  a  matter  of 

fact,  such  a  hard  :  oivtrun.t  may  ;  .  :  ’  an  i  nnp.it.  ihie 

constraint  set,  as  demons  t  rated  ry  •  joi  lowing  example 
from  the  spectral  estimation  arena. 

Suppose  a  small  number  of  low-lag  values  or  <•»  covariance 
sequence  is  available,  and  extrapolation  ■  l  t  hew  i  .  ues  i  red 
in  order  to  increase  spectral  resolution.  A  thooret  real 
constraint  on  the  spectral  density  function  would  be  nor.no- 
gativity.  Now  suppose  that  due  to  measurement  noise  the 
zero  lag  covariance  value  is  not  the  largest  in  magnitude. 

This  now  violates  one  of  the  necessary  properties  for 
covariance  sequences,  and  results  in  an  incompatibility  of 

the  given  covariance  sequence  segment  with  any  spectral 
density  function. 
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In  order  to  accomodate  errors  :  r.  the  :  1  ab.e 

measurements,  we  introduce  the  following  so:;  .  or.,  i  r  aim  a:, 
frequency  domain  measurement  constraint. 


if  I  I  F-Y  I  iM  -  6  -  f- 


cff 


Y  +  |  (F-Y)  over  M 


over  M 


(20) 


if  I .F-Y| ,M 


Herein,  t  indicates  the  tolerable,  and  5“  indicates  the  ac¬ 
tual,  mean  square  difference  between  reconstruction  ana  mea¬ 
surements.  These  evaluations  are  made  over  the  ir.eu.,si  *-ment. 
domain  M  only.  If  corrective  action  takes  place,  the  fol¬ 
lowing  relationship  holds 


1 ICfF-Y| 

•m 

"  MY  ♦  f  (F-Y)-Yi  lM 

= 

1  MF-Y||M  -  E 

Therefore , 

we 

1  inu , 

1 ICfF-Y| 

i2 

=  MF-YM^c  *  |'l\,F-Y||-; 

- 

If  *  ~ Y  i  '  Mc  +  7*  1  I  ^  '  :  M 

a 

1  1 F-Y |  |  2 

which  says 

that  the  iterate  either 

remains  unchar 

moves  closer 

to  the  measurements. 

Our 

soft  mea 

constraint  leaves  well  enough  alone,  i.e.  if  the  iterative 
process  comes  clo. e  enough  to  Y,  with  respect  to  the  expect- 
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ed  noise  level,  then  the  1  torut  ••  i.  cons: do i  ei  uc>.  vpv, . 
An  alternative  soft  constraint  would  be  one  that  placet  a 
tolerated  maximum  on  deviations  from  the  measurements,  for 
each  measurement  sample  individually.  Note  that  in  any  case 
a  weighted  norm  is  easily  accomodated,  so  that  known  ampli- 
tude-and/or  frequency-dependent  noise  j  1.1  ■.<:  mti  t .  •  r.  may  re  .d;  - 
ly  be  incorporated. 

Several  constraints  apply  in  the  scene  domain.  As  the 
solution  to  our  problem  is  a  picture,  compatible  with  our 
information,  a  nennegat 1 v : ty  operator  up pi ies ,  which  is  net 
softened.  The  scene  is  also  of  limited  extent  and  it  is  ar¬ 
gued  that  softening  may  be  beneficial  here .  A  border  region 
B  is  proposed  (similarities  with  I'?1),  to  serve  as  a  transi¬ 
tion  region  between  the  known  extent  K,  and  the  zeroed  re¬ 
gion  Z.  The  soft  scene  limitation  operator  S  then  becomes 

over  K 

over  B  (!.':) 

over  Z 

where  0  n  gD  <  1  is  a  function  defined  ov**r  the  horde:  re- 

gion,  which  provides  a  smooth  transition  to  zero. 

In  the  border  region  itself  one  could  elect  not:  to  make 
any  substitutions,  and  then  monitor  the  energy  in  the  border 
region.  This  information  can  potentially  be  used  to  accel¬ 
erate  the  convergence  of  the  algorithm.  We  hope  to  investi¬ 
gate  this  issue  in  the  subsequent  stage  of  research. 
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C.  Convergence  Issues 

In  order  to  assess  the  prospect:;  :  c  j  c  or.'.’er  g-no--  .  ;  (  ]  » 

the  contraction  properties  of  several  operators  must  be 
evaluated  [6,81.  For  the  scene  limitation  operator  -S,  via 
Parseval's  relationship 


SF_1E'1-SF"1Fj 


|f,-f  .  |“dt  ’  — 


_  i  - 1 

F  F  -  F  1  F 


:  - 1 


3  ' 

(13  ) 


whe  re 


U  I  f ,  - f ,  |  “dt  +  L  ( 1  -  I  gB  i  " .)  I  t  ,  -  i 


2  ,  2"  i2'*1  ^ 

r  =  1 - 


>1!  JB1 


I  I  F  1Fi  -  r_1F  | 


(24) 


The 


expression  for  r“  is  nonnegative  and  smaller  or  equal 


one,  where  equality  occurs  if 


union  of  B  and 


and  f  are  euua 1  over 


t;.e 


The  operator  5  is  there! ore  nonexpans ive . 
The  nonnegativity  operator  can  similarly  be  shown  to  be 
nonexpansive . 


I  I  PF " 1 - TF _  1  F 


,  I  IjPtd-fY  i“  lt 


if  '  "dt. 

i  3 


|  F~  1 F  -  F  1 F 
i  j 


(2r-) 


Equality  holds  if  the  functions  sign  (  t  )  and  sign  1 1.  )  ate 
equal  almost  everywhere  1  -  i  ,  n  |.  The  composite  scene  domain 
operator 


4 


-Li- 


I 


Ce  =  PS 


=  SP 


is,  consequently,  a  nonexpansive  operator. 

The  soft  frequency  domain  measurement,  constraint  yields, 


I  C  -  F  -C-F . 

r  f  i  c  j 


||  r  F  -C  c  i  *  -  1  r  t.- 

1  1  c  F  i  F  I  '  1  v  F  ,  F 


<  "Fi-FJMK 


^  I  IF  -F; 

1  j 


The  inequality  arises  because  the  constraint  leaves  the  com¬ 
ponents  of  F^  and  fF  in  M  alone,  or  moves  cue  or  both  compo¬ 
nents  toward  Y,  and  consequently  closer  together.  The  com- 
c 

ponents  in  M  are  always  unchanged.  We  find  the  measurement 
operator  C^,  therefore,  to  be  a  nonexpansive  operator  also. 
Rests  us  to  evaluate  the  operator  (I-uT). 


Fi-uTFi-F]  +'j7Fj  I  > 


I  f  -uF .  -  f  .  +uF  i  i ;;  <  j  |  V  -F 


till 


1  J  MC 


| |Fi-Fj I |  for  0  ' 


Equality  holds  if  F^  equals  F  over  M. 

The  conditions  for  equality  in  (24),  (2r>),  (29),  and  (32) 

could  all  be  met  at  once,  and  consequently  we  find  (I-yT)C 


to  be  a  non- expansive  operator.  As  a  result  (IB)  may  have 
many  fixed  points.  It  should  be  recognised  tna  t  i  n  m  ew  oi 
the  noisy  measurements,  a  unique  solution  cannot  be  expect¬ 
ed.  Any  solution,  however,  must  satisfy  all  the  constraints. 
The  size  of  the  solution  set  can  possibly  be  reduced  by 
decreasing  c  in  (20) ,  but  his  moves  the  problem  towards 
incompatibility  of  the  available  information. 

D.  Algorithm  Simplification 

Using  the  soft  frequency  domain  constraint  m  (20),  and 
the  composite  scene  domain  cc  i.sti  aim  m  (20),  in  tir.-  algor¬ 
ithm  formulation  of  (18),  yields  the  following  algorithm. 


Fo  =  '-Y 

Fk.i  =  “Y  *  !H-fii ysrh-K  on 

This  algorithm  is  represented  m  Figure 
Let's  define 


CFF^SF  Fk 


(  34) 


so  that  the  right  hand  side  of  (33)  c  an  be  written  and  sim¬ 
plified  as  follows. 


Fr  -  UT F y  +  yY  -  F  over  Mc 


(3c) 


=  (I-u)Fr  -i  uY  over  M 


(36) 


The  latter  equality  can  be  rewritten  as 
Fk-yTFk+pY  -  Y+(l-y)(Fk-Y) 


(37) 
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If  we  now  substitute 


t 

6 


(37  ) 


for  (  1  -  u  )  t  i. 


seen  to  represent  the  sort  frequency 
(20).  As  ~  is  smaller  or  equal  one 
choice  for  p : 


..3  1  and 
.  r  ana  n.  con:.  :.r  amt 
the  corresponding 


U  1  "  6  (38) 

is  seen  to  satisfy  the  convergence  condition:.,  ox  i  3.':  )  .  A.-:  a 

result  of  making  this  particular  choice  :<•»  w<-  wru.u  nave 

m  Figure  2,  two  identical  rcnrecux.'  u-tavy  ;u 

constraint  operators.  The  algor:  t.h:n  thuretote  simplifies 
because  the  dot-enclosed  algorithm  part  now  has  the  effect 
of  an  identity  operator,  due  to  the  equality 


C  • 

F 


C 


F 


1  -c 


(  3  d  ) 


The  resulting  algorithm  simply  consists  oi  transforma¬ 
tions  from  frequency  domain  to  :.eeno  domain  and  vi  ue  versa, 
with  an  application  of  the  respective  constraint.  operators. 

Note  that  if  our  tolerable  noire  level  equals  uero, 
c  =  0  in  (20),  then  (20)  as  well  as  (38)  cause  the  algorithm 
to  degenerate  to  the  traditional  hard-constraint  algorithm, 
in  which  the  iterate  is  replaced  with  the  measurements  over 
the  measurement  domain. 


( 
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4 .  SIMULATIONS 

For  each  of  the  blocks  in  Figure  2  commands  were  written 
in  RATFOR,  for  subsequent  running  on  the  VAX  ll/78o  m 
VFI&SU's  SDA  lab.  A  command  is  like  a  self-contained 
subroutine  called  from  the  keyboard,  and  usually  operating 
on  existing  files.  Each  iteration  of  the  algorithm 
therefore  takes  four  commands.  One  more  command  was  written 
to  take  care  of  all  initializations.  Codes  for  these 
commands  can  be  found  in  the  Appendix. 

Our  scene  domain  consists  of  a  64  x  64  pixel  array. 
After  generating  a  checkerboard  pattern  the  scene  domain  was 
constrained  to  the  center  25  x  20  (hor  x  vert)  pixels,  which 
therefore  constitute  the  region  of  known  extent  K.  The 
resulting  scene  limited  picture  was  then  transformed  and  the 


measurement 

domain  M 

was  chosen  to  be¬ 

the  15  x  19 

(hor  x 

vert)  area 

centered 

at  zero  frequency 

/ .  The  mean 

squa re 

value  of  the  frequency  domain  signal  turned  out  to  be  27  x 
10^.  Noise  can  be  added  to  these  frequency  domain 
measurement  values.  We  used  a  noise  variance  of  40  x  103  . 
The  scene  domain  constraint  implemented  in  these  simulations 
left  a  scene  of  27  x  22  (h-.-r  x  vert)  pixels,  so  that  the 
border  region  B  is  1  pixel  wide  on  each  side. 

At  present  the  set  of  commands  has  to  be  entered 
sequentially  in  order  to  effect  a  number  of  iterations  of 
the  algorithm.  This  time  consuming  process  can  be 
alleviated  in  the  future  by  rewriting  the?  commands  into  a 


■  i 

i 


-18- 


single  command,  which  will  ’.hen  l,i  liitute  a  ir ok- 

comprehensive  study  or  the  anvenj-no;  behavior  i  or  re¬ 
numbers  of  iterations  of  the  pi  euent  algorithm.. 

Consequently  our  present  numerical  results  are  limited  to  7 
iterations  of  the  algorithm.  Tin..  i  nmted  experiment 

however,  demonstrates  what  ,5  expected  to  bo-  r  he 
characteristic  behavior  of  the  algorithm. 

Starting  from  the  initial  condition  without  noise  and 

going  through  7  iterations  will  not  lean  to  any  uramatic 

results.  We  illustrate,  never  the  It-.  s ,  the  operational 

condition  of  the  algorithm  m  Figure  J.  where  a  cur;  tar.tial 

improvement  ir.  resolution  is  ol-t  n :  n*-d  alter  on  1  y  ..even 

iterations,  when  compared  tu  a  ?  1  mt.>  1  •.•  in  /er  tram  ic::r,  ol 

the  measurements.  To  illustrate  the  ...  i.<u  ten  sti  c  behav;  .  r 

.  :  tiie  soft-con.-tr.i :  nt  a  1  g-.-r  1  t  hm  t  i-,;:  :  :  r.  •  y  7 

iterations,  the  next  ext  on  men t  ,  w  .  t  h  :.  ■  : ; y  rnu.-a./ar  eni>-:it.;. 

(CNR  --  28  dS)  starts  trom  the  o*.  n-*i  --n:.  her  the  true 

solution  be  the  initial  iterate,  which  should  therefore 

constitute  a  feasible  solution  foi  *  n- •  a ;oori thm.  The 

iterate  F,  in  the  traditional  aim  •  :  t  nr  with  the  hard 

k 

frequency  domain  constraint  then  move* ..way  :  rom  the  true 
solution  F,  thereby  indicating  that  th. :  algorithm  cannot 
possibly  converge  to  the  true  r.olut  ion.  Tin oftei  t  w«; 
actually  ooservable  on  the  TV  monitor,  al  though  it  doer-  no 
show  in  the  photograph  of  Figure  4.  W»-  therefore  derived 
numerical  indicators  of  the  behavior  of  the  reconstructions, 
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Upper  Left;  Original  Scene 

Upper  Right;  Inverse  of  Measurements  without  Noise 

Lower  Left;  Iterate  7  before  Scene  Domain  Constraint  Operator 

Lower  Right;  Iterate  7  after  Scene  Domain  Constraint  Operator 

Figure  3.  Algorithm  Operation  for  Noiseless  Measurements, 

with  Inverse  of  Measurements  as  Initial  Interate. 


Upper  Left:  Original  Scene 

Upper  Right:  Inverse  of  Measurements  with  Noise 

Lower  Left:  Iterate  ?  before  Scene  Domain  Constraint  Operator 

Lower  Right:  Iterate  7  after  Scene  Domain  Constraint  Operator 

Figure  U.  Algorithm  Operation  for  Noisy  Measurements, 
with  Original  Scene  as  Initial  Interate. 


I 

and  these  RMS  values  are  represented  :n  Fiuure  z> .  T:. 
iterate  in  the  soft  constraint  „*••?  •!  :  thrr.  with  the  no  i . 
variance  underestimated  at  10, COG  ( the  raise  variance  evua ! 
40,000)  is  seen  to  also  move  away  from  the  true  solution 
but  much  slower.  Finally,  it  is  clear  that  if  a  goo 
estimate  of  the  noise  variance  is  ava.iabie,  then  th 
iterate  will  not  move  away  from  the  ti  ue  solution,  Lecaus 
they  are  close  enough  already  to  satisfy  the  sol 
constraint.  The  specific  effects  of  estimated  and  actua 
noise  variance  will  be  investigated  .  r.  .«  ...ubsequent  lesears 
effort. 


_  )  >  _ 
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# — FFTDRV  DRIVER 

« 

♦IDENTIFICATION 

#  TITLE  FFTDRV 

#  AUTHOR  AA(LOUIS)  3EEX 

«  VERSION  A. 01 

#  DATE  15  SEPT  1982 

#  LANGUAGE  RATFOR 

#  SYSTEM  VAX-11 

# 

♦PURPOSE 

♦  DRIVER  FOR  COMMAND  THAT  PERFORMS  THE  RADIX-2  FAST 

FOURIER  #  TRANSFORM  OF  A  REAL  IMAGE 

♦ 

♦ENTRY  POINT 

♦  FFTDRV (WORK ,ERRET) 

♦ 

♦ARGUMENT  LISTING 

#  WORK  I NT  WORK  ARRAY 

#  ERRET  INT  ALTERNATE  ERROR  RETURN 

# 

♦INCLUDE  FILES/COMMONS 

♦  MACA1  INCLUDE  GIPSY  TOKEN  DEFINITIONS  FOR  A1 

CHARACTER  #  GIPCOM  INCLUDE 

♦  ERRET  INCLUDE  INCLUDE  FILE  FOR  COMMON  ERROR 

♦ 

♦ROUTINES  CALLED 

#  PPUSH  PUSHES  PROGRAM  NAME  INTO  ERROR  STACK 

(GIPSY)  ♦  PPOP  POPS  PROGRAM  NAME  FROM  ERROR  STACK 

(GIPSY)  ♦  RDKINL  OPEN  A  SIF  AND  FILL  IN  THE  IDENT 

BLOCK  #  CLOSE  CLOSES  A  FILE (GIPSY) 

♦  FFTALG  COMPUTES  THE  (INVERSE)  FAST  FOURIER 

TRANSFORM  (USER)  ♦ 


****************************************************************** 

# 

INCLUDE  MACA1 
SUBROUTINE  FFTDRV ( WORK, *) 

IMPLICIT  INTEGER  (A-Z) 

INCLUDE  GIPCOM 

INCLUDE  ERROR 

INTEGER  IDENT ( . IDLENGTH ) 

DIMENSION  WORK (.ARB)  ,  , 


EQUIVALENCE (NPPL, IDENT ( . IDNPPL) ) 

EQUIVALENCE (NLIN, IDENT ( . IDNLINS) ) 

EQUIVALENCE (NCOL, IDENT (. IDNCOLS) ) 

EQUIVALENCE (NROW, IDENT (. IDNROWS) ) 

EQUIVALENCE (NBND, IDENT ( . IDNBNDS) ) 

EQUIVALENCE (MODE, IDENT ( . IDMODE) ) 

# 

CALL  PPUSH ( ' FFTDRV ' ) 

# 

#  OPEN  INPUT  FILE 

# 

CALL  RDKINL ( FDI1 , IDENT, .OLD, IEV,%9999) 

CALL  CLOSE (FDI1) 

# 

#  CHECK  INPUT  FILE 

* 

IF (MODE~= . REALMODE)  GOTO  9000 
M=1 

NLIN2=NLIN/2 
WHILE (NLIN2>1) 

$( 

M=M+1 

NLIN2=NLIN2/2 

$) 

IF(NLIN,/'-2**M)  GOTO  9020 
N=1 

NPPL2=NPPL/ 2 
WHILE (NPPL2>1) 

$( 

N=N+1 

NPPL2=NPPL 2/2 
$) 

IF  ( NPPL”'  =  2**N)  GOTO  9020 
NMAX=MAX (NPPL, NLIN) 

# 


# 

NXT=1 

ILIN1=GETWP (NXT, . REALMODE , NPPL) 

ILIN2=GETWP(NXT, .REALMODE, NPPL) 

XAR=GETWP ( NXT , . REALMODE , 2  *NMAX ) 

FTAR=GETWP (NXT, . REALMODE, 2*NLIN*NPPL) 

* 

IF ( . OK~=OSALOC ( NXT) )  GOTO  9010 

* 

CALL  COMTIN (%9999) 

* 

#  CALL  FFT  ROUTINE 

# 

CALL  FFTALG (FDI1,FD01, NBND, WORK (IL INI) ,WORK(ILIN2) , WORK (XAR) , 
WORK ( FTAR) , NL IN, NPPL , NMAX , IEV,%9999) 

* 

CALL  PPOP 

* 

RETURN 

* 

#  ERROR  CONDITIONS 

9000  CONTINUE 

I 


I 


#  ILLEGAL  I 'AT  A 

IEV=-2012 
GOTO  9999 

# 

9010  CONTINUE 
IEV=OSGIEV ( IEV) 

GOTO  9999 

# 

#  ILLEGAL  ARRAY 

9020  CONTINUE 
IEV=-5004 
GOTO  9999 

# 

9999  CONTINUE 
CALL  CLOSE (FDI1) 

CALL  CLOSE (FDOl) 

# 

RETURN  1 
END 


MODI'; 


SIZE 


A-  I 


I/O  t,  ALGORITHM 


# — FFTALG 
# 

# IDENTIFICATION 

#  TITLE 

#  AUTHOR 

#  VERSION 

#  DATE 

#  LANGUAGE 

#  SYSTEM 

# 

♦PURPOSE 

#  TO  COMPUTE  THE 

#  OF  A  SIF  IMAGE 

# 

♦ENTRY  POINT 


FFTALG 

AA(LOUIS)  BE EX 
A. 01 

15  SEPT  1982 

RATFOR 

VAX-11 


FAST  FOURIER  TRANSFORM 
(1  BND  IF  REAL,  2  ENDS 


(RADIX-2) 

IF  COMPLEX) 


FFTALG (FDI,FDO, NBND, ILIN1 , XL IN2 , XAR , FTAR , ML I N, NPPL , UMAX , IEV, * ) 
♦ARGUMENT  LISTING 


♦ 

FDI 

INPUT  FILE  DESCRIPTOR 

♦ 

FDO 

OUTPUT  FILE  DESCRIPTOR 

♦ 

NBND 

NO  OF  BANDS  IN  IMAGE  TO  BE  PROCESSED 

# 

ILIN1 

LINE  BUFFER  FOR  REAL  PART  OF  INPUT  IMAGE 

♦ 

ILIN2 

LINE  BUFFER  FOR  IMAGINARY  PAR'!’  OF  IMAGE 

♦ 

XAR 

WORK ARRAY 

♦ 

FTAR 

WORK ARRAY 

♦ 

NLIN 

NO  OF  ROWS  IN  THE  IMAGE 

# 

NPPL 

NO  OF  COLUMNS  IN  THE  IMAGE 

♦ 

NMAX 

MAX  OF  NPPL  AND  NLIN 

♦ 

INV 

INVERSE  TRANSFORM  IF  TRUK 

♦ 

IEV 

INTEGER  EVENT  VARIABLE 

# 

# 

♦INCLUDE 

ERRET 

ALTERNATE  RETURN 

:  FILES/COMMONS 

♦ 

U 

MACAl 

INCLUDE  (GIPSY) 

If 

♦ROUTINES  CALLED 

♦ 

PPUSH 

♦ 

PPOP 

♦ 

CLOSE 

♦ 

CPYIDK 

# 

COPYDS 

♦ 

RREAD 

♦ 

RWRITE 

♦ 

RX2FFT 

♦ 

DSCNAM 

♦ 


★  ★★a********************'*******************-********'************ 

# 

INCLUDE  MACAl 
SUBROUTINE 

FFTALG (FDI1 , FDOl , NBND , IL I Nl , IL IN2 , XAR , FTAR , ML I N , NPPL , NMAX , T  EV, * 
IMPLICIT  INTEGER  (A-Z)  CHARACTER  FDI 1  (  .  FDLENGTI1 )  , 

FDOl ( .FDLENGTH)  INTEGER  IDENT ( . IDLENGTH) ,  JDENT ( . IDLENGTH) 

REAL  ILIN1 (NPPL) ,  ILIN2(NPPL),  NORM 
COMPLEX  XAR ( NMAX ) ,  FTAR ( NLIN, NPPL) 

COMPLEX  Z 
LOGICAL  INV 

# 

INCLUDE  TTCOM 

♦ 


•  *4 


CALL  PPUSH ( ' FFTALG 1 ) 

# 

NMAX=MAX (NLIN, NPPL) 

F0RWARD=1 

INVERSE=-1 

DFLT=1 

CALL  RNGETI ( 'FORWARD (1)  OR  INVERSE ( - 1 ) 

TRANSFORM?. ' , INVERSE , FORWARD , DFLT ,  VAR, IEV,%9000) 

INV=.TRUE. 

IF ( VAR==1 )  INV=. FALSE. 

#  OPEN  INPUTF1 LE 

CALL  CPYIDR(FDIl , IDENT, . OPNTMP, IEV,%9000) 

# 

#  WRITE  DESCRIPTOR  RECORD  TO  TEf'l 

FILE  CALL  DSCNAM ( 1 FFTALG ' , IEV, % 9000) 

# 

#  SET  UP  OUTPUT  FILE  PARAMETERS 
JDENT ( . IDNPPL) =IDENT( . IDNPPL) 

JDENT ( . IDNLINS) =IDENT( . IDNLINS) 

JDENT ( . IDNCOLS) =IDENT { . IDNCOLS) 

JDENT ( . IDNROWS) = IDENT ( . IDHROWS) 

JDENT ( . IDNBNDS) =2 
JDENT ( . IDMODE) =2 

# 

IAN  INVERSE  TRANSFORM  MUST  YIELD  A  REAL  IMAGE  FOR  OUR 
APPLICATION:  #NO.  OF  OUTPUT  BANDS  EQUALS  1 
NOB=2 

IF(INV)  NOB=l 
JDENT ( . IDNBNDS) =NOB 

# 

#  OPEN  OUTPUT  FIEF 

CALL  C0PYDS(FD01 , JDENT, IEV, %  9000) 

# 

DO  BLKN=1 , NLIN 
$( 

CALL  RREAD (FDI1,ILIN1,1,BLKN, IDENT, .WAIT, 1EV,%9000) 

IF (NBND==2 ) 

$( 

CALL  RREAD  (FDI1,ILIN2,2, BLKN,  IDENT,  .WAIT,  1EV, *9000) 

$) 

ELSE 
$  ( 

DO  1=1, NPPL 
$( 

ILIN2 (I) =0. 

$) 

$) 

IF(INV) 

$( 

DO  J=1 , NPPL 
$( 

XAR ( J ) =CMPLX ( IL INI ( J ) , IL IN2 ( J ) ) 

XAR ( J ) =CONJG ( XAR ( J ) ) 


$( 

DO  J=1 , NPPL 
$  ( 

XAR(J) =CMPLX ( ILIN1 ( J) ,ILIN2(J) ) 

$) 

$) 

CALL  RX2FFT (XAR, NPPL) 

DO  J  =  1 , NPPL 

$< 

FTAR(BLKN, J) =XAR(J) 

$) 

§) 

NORM=FLOAT ( NLIN*NPPL ) 

DO  J  =  1 , NPPL 
$( 

DO  1=1 , NL IN 
$( 

XAR ( I )  =  ( FTAR ( I ,  J )  ) 

$) 

CALL  RX2FFT (XAR, NLIN) 

IF(INV) 

$( 

DO  1=1, NLIN 
?  ( 

FTAR ( I , J) =CONJG ( XAR ( I ) ) / NORM 
$) 

$) 

ELSE 
$  ( 

DO  1=1, NLIN 
$  ( 

FTAR ( I , J ) =X AR ( I ) 

$) 

S) 

$) 

# 

DO  BLKN=1 , NLIN 
$  ( 

DO  1=1, NPPL 
$  ( 

Z=FTAR ( BLKN, I ) 

ILIN1 ( I ) =REAL ( Z) 

ILIN2 ( I ) =AIMAG ( Z) 

$) 

CALL  RWR ITE ( FDOl , I L INI , 1 , BLKN , JDENT , . WA I T , I EV , %  9 0 0  0 ) 
IF ( NOB==l )  GOTO  1000 

CALL  RWRITE (FDOl , ILIN2 , 2 , BLKN, JDENT, . WA IT , I EV , % 9 0 0 0 ) 
1000  CONTINUE 
$) 

CALL  PPOP 
CALL  CLOSE (FDI1) 

CALL  CLOSE (FDOl) 

RETURN 

# 

# 

#  ERROR  CONDITIONS 

9000  CONTINUE 
CALL  CLOSE (FDI1) 


CALL  CLOSE (FDOl) 

* 

9999  CONTINUE 
9090  RETURN  1 
END 


# — RX2FFT  THIS  IS  A  SUBROUTINE  CALLED  IN  PFTALG.RAT 

SUBROUTINE  RX2FFT(X,N) 

COMPLEX  X (N) ,  U,  W,  T 

NV2=N/2 

NM1=N-1 

M=1 

WHILE (NV2>1) 

$( 

M=M  +  1 
NV2=NV2/2 
?) 

NV2=N/2 
J  =  1 

DO  1=1, NM1 
$  ( 

IF ( I >=J)  GOTO  5 
T=X (J) 

X(J)=X(I) 

X ( I ) =T 

5  K=NV2 

6  IF ( K >  =  J )  GOTO  7 
J=J-K 

K=K/2 
GOTO  6 

7  J  =J +K 
$) 

PI=4 . *ATAN ( 1 . ) 

DO  L  =  1 , M 
$( 

LE=2**L 
LEI =LE/ 2 
U=CMPLX(1. ,0.  ) 

FLE1=FL0AT(LE1) 

CARG=COS ( PI/FLE1 ) 

SARG=-1.*SIN(PI/FLE1) 

W=CMPLX (CARG , SARG ) 

DO  J=1 , LEI 
$  ( 

DO  I=J,N,LE 
$  ( 

IP= I+LE1 
T=X (IP) *U 
X ( IP) =X ( I) -T 
X ( I ) =X ( I )  +T 
$) 

u=u*w 

$) 

$> 

RETURN 

END 


This  subroutine  is  an  adaptation 
algorithm  for  decimation-in-time 


of  tlie  Cooley,  Lewis,  and  Welch 
radix-2,  in-place  FFT. 


#  — SCNCTD  DRIVER 

♦ 

# IDENTIFICATION  BEEX ,  RATFOR,  VAX-11 

* 

# PURPOSE 

#  DRIVER  FOR  COMMAND  THAT  ZEROES  A  SCENE 

#  OUTSIDE  OF  A  SPECIFIED  WINDOW 

# 

♦ENTRY  POINT  SCNCTD (WORK , ERRET) 

# 

♦INCLUDE  FILES/COMMONS 

♦  MACA1 

♦  GIPCOM 

♦  ERRET 

♦ 

♦ROUTINES  CALLED 

♦  PPUSH 

♦  PPOP 

♦  RDKINL 

♦  CLOSE 

♦  SCNCTA  ALGORITHM  SECTION  OF  SCNCTT  COMMAND 

♦ 

********************************  ********************************* 
♦ 

INCLUDE  MACA1 
SUBROUTINE  SCNCTD (WORK ,* ) 

IMPLICIT  INTEGER  (A-Z) 

INCLUDE  GIPCOM 

INCLUDE  ERROR 

INTEGER  IDENT( . IDLENGTH) 

DIMENSION  WORK (.ARB) 

♦ 

EQUIVALENCE  (NPPL, IDENT ( . 1DNPPL) ) 

EQUIVALENCE  (NLIN, IDENT ( . IDNLINS) ) 

EQUIVALENCE  (NCOL, IDENT ( . IDNCOLS) ) 

EQUIVALENCE  (NROW, IDENT ( . IDNROWS) ) 

EQUIVALENCE  (NBND, IDENT ( . IDNBNDS) ) 

EQUIVALENCE  ( MODE , IDENT ( . IDMODE) ) 

♦ 

CALL  PPUSH ( ' SCNCTD ' ) 

I 

♦  OPEN  INPUTF ) LE 

♦ 

CALL  RDKINL (FDI1 , IDENT, .OLD, IEV,%9999) 

CALL  CLOSE (FDI1) 

♦ 

♦  •  CHECK  INPUTFILE 

♦ 

IF (MODE~= . REALMODE)  GOTO  9000 

♦ 

NP=NPPL*NLIN 

N=1 

NP2=NP/2 
WHILE (NP2>1) 

$  ( 

N=N+I 

NP2=NP2/2 

$) 

IF(NP~=2**N)  GOTO  9020 
IF ( NBND*=1 )  GOTO  9020 

#  ,  A-«. 

NXT=1 


ILIN1=GETWP (NXT, . REALMODE, NPPL) 

IF ( . OK~=OSALOC ( NXT) )  GOTO  9010 

# 

CALL  COMTIN (%9999) 

# 

#  CALL  ALGORITHM  SECTION 

# 

CALL  SCNCTA ( FDI1 , FDOl , WORK ( I L INI ) ,NPPL,NLIN, JEV,%9999) 

# 

CALL  PPOP 

# 

RETURN 

#  ERROR  CONDITIONS 

9000  CONTINUE 

#  ILLEGAL  DATA  MODE 
IEV=-201 2 

GOTO  9999 

# 

9010  CONTINUE 
IEV=OSGIEV ( IEV) 

GOTO  9999 

#  ILLEGAL  ARRAY  SIZE 
9020  CONTINUE 

IEV=-5004 
GOTO  9999 

# 

9999  CONTINUE 
CALL  CLOSE (FDI1) 

CALL  CLOSE (FDOl) 

# 

RETURN 1 
END 


i 

I 


I 


# — SCNCTA  ALGORITHM 

# 

# IDENTIFICATION  BE EX ,  RATE OR,  VAX-11 

# 

IPURPOSE 

#  ALGORITHM  SECTION  FOR  COMMAND  THAT  ZEROES 

#  A  SCENE  OUTSIDE  OF  A  SPECIFIED  WINDOW 

# 

#ENTRY  POINT 

#  SCNCTA ( FDI ,FDO, ILIN1 , NPPL, NLIN , IEV, *) 

# 

# ARGUMENT  LISTING 

#  FDI  INPUTFILE  DESCRIPTOR 

#  FDO  OUTPUTFILE  DESCRIPTOR 

#  ILIN1  LINE  BUFFER 

#  NPPL  NO  POINT  PER  LINE 

#  NLIN  NO  LINES 

#  IEV 

# 

# 

INCLUDE  MACAl 

SUBROUTINE  SCNCTA  ( FD 1 1 ,  FDOl  ,  ILIN.1 ,  NPPL, NLIN,  IEV,  *) 

IMPLICIT  INTEGER  (A-Z) 

CHARACTER  FD I 1 ( . FDLENGTH ) ,  FDOl ( . FDLENGTH ) 

INTEGER  IDENT ( . IDLENGTH) ,  JDENT ( . 1 DLENGTH ) 

REAL  ILIN1 ( NPPL) 

* 

INCLUDE  TTCOM 

# 

CALL  PPUSH( 'SCNCTA' ) 

# 

HLOW=2 
HHI=NPPL-1 
DFLT=NPPL/ 10 

CALL  RNGETI ( 'HORIZONTAL  WINDOW 
SIZE?.  '  ,  BLOW,  HH I  ,DFLT,  WXI.EN ,  I EV,  I.  ROOD  )  VLO=2 
VHI=NL IN-1 
DFLT=NLIN/10 

CAl.L  RNGETI  ( 'VERTICAL  WINDOW 

SIZE?.  '  , VLO, VHI,DFLT, WYLEN,  IEV,%9000)  IF  (WXLEN~=  ( WXLKLV  2)  *2) 
WXLEN=WXLEN+1  IF ( WYLEN (WYLEN/ 2 ) *2)  W Y  L  H  N = W Y  L  E N  + 1 


I  OPEN  INPUTFILE 

i 

CALL  CPYIDR ( FDIl , IDENT, . OPNTMP, IEV,%9000) 
CALL  DSCNAM ( ' SCNCTA' ,IEV,%9000) 


DO  1=6,19 
$( 

IF  (I  .EQ.  6  .OR.  I  .EQ.  7  .OR.  I  .EQ.  13  .OR.  I  . EQ.  14  .OR. 
I  .EQ.  17  .OR.  I  .EQ.  19)  $( 

X= IDENT ( I ) 

JDENT ( I ) =X 

$) 

$) 

#  OPEN  OUTPUTFILE 
CALL  C0PYDS(FD01, JDENT, IEV,%9000) 

I 

IF (WXLEN<  =0/WXLEN>=NPPL/WYLEN<  =  0/WYLEN>  =  NLIN)  GOTO  9000 

# 

XWM1= (NPPL-WXLEN) /2 


4*  4*  4* 


XWB=XWMl->  1 

XWEP1=XWB+WXLEN 

XWE=XWEP1-1 

# 

YWM1= ( NLIN-WYLEM) /2 
YWB=YWM1+1 
YWEP1=YWB+WYLEN 
YWE=YWEP1 -1 

# 

DO  I=1,NPPL 
$  ( 

ILIN1 ( I ) =0 . 

$) 

DO  BLKN=1 ,  YWM1 
$  ( 

CALL  RWRITE (FDOl , ILIN1 , 1 , BLKN, JDENT, . WAIT, IEV, % 9000) 

$) 

DO  BLKN=YWEP1 , NLIN 
$  ( 

CALL  RWRITE  (FDOl,  ILIN1 , 1 ,  BLKN ,  JDENT,  .WAIT,  IFV,  7.9000) 

$) 

DO  BLKN=YWB, YWE 
$  ( 

CALL  RREAD (FDIl , ILIN1 , 1 , BLKN, IDE NT, .WAIT, IEV, t 9000) 
DO  1=1 ,XWMl 
$  ( 

ILIN1 ( I ) =0 . 

S) 

DO  I=XWEP1 , NPPL 
$  ( 

IL INI ( I) =0  . 

5) 

#  POSITIVITY  CONSTRAINT 
DO  I=XWB , XWE 

$  ( 

IF ( ILIN1 ( I ) <  0  . )  ILIN1 (I) =0. 

S) 

CALL  RWRITE  (  FDOl ,  JLINl,!  ,  BLKN,  JDENT,  .WAIT,  IEV, V, 9000) 

$) 

# 

CALL  PPOP 
CALL  CLOSE (FDIl) 

CALL  CLOSE (FDOl) 

RETURN 


ERROR  CONDITIONS 


9000  CONTINUE 
CALL  CLOSE ( FDIl ) 
CALL  CLOSE (FDOl) 

# 

9999  CONTINUE 
9090  RETURN  1 
END 


#  --FRQCTD  DRIVER 

# 

# IDENTIFICATION  BEEX ,  RATFOR,  VAX-11 

« 

#  PURPOSE 

#  DRIVER  SECTION  FOR  COMMAND  THAT  IMPOSES  A 

#  FREQUENCY  DOMAIN  MEASUREMENT  CONSTRAINT 

# 

# ENTRY  POINT  FRQCTD (WORK, ERRET) 

# 

# INCLUDE  FILES/COMMONS 

#  MACA1 

#  GIPCOM 

#  ERRET 

# 

# ROUTINES  CALLED 

#  PPUSH 

#  PPOP 

#  RDKINL 

#  CLOSE 

#  FRQCTA  ALGORITHM  SECTION  OF  l'RQCTC  COMMAND 

******************************************************* 

# 

INCLUDE  MACA1 
SUBROUTINE  FRQCTD (WORK ,* ) 

IMPLICIT  INTEGERCA  -  Z) 

INCLUDE  GIPCOM 
INCLUDE  ERROR 

INTEGER  IDENT (. IDLENGTH) ,  LDENT  (  .  IDLENGT1I ) 

DIMENSION  WORK (.ARB) 

# 

EQUIVALENCE  (NPPL,  IDENT  (  .  IDNPPI.)  ) 

EQUIVALENCE  (NLIN, IDENT ( . 1DNL1NS) ) 

EQUIVALENCE  (NCOL, IDENT ( . IDNCOLS) ) 

EQUIVALENCE  (NROW, IDENT (. 1DNROWS) ) 

EQUIVALENCE  ( NBND, IDENT ( . IDNBNDS) ) 

EQUIVALENCE  ( MODE , IDENT ( . IDMODE) ) 

# 

CALL  PPUSH ( 'FRQCTD' ) 

#  OPEN  IMPUTE  ILLS 
CALL  RDKINL (FDI1, IDENT, .OLD, IEV,%9999> 

CALL  CLOSE (FDI1) 

CALL  RDKINL (FDI2, LDENT, .OLD, IEV,%9999) 

CALL  CLOSE (FDI2) 

#  CHECK  INPUTFILKS 
IF (MODE~=. REALMODE)  GOTO  9000 

#NEED  TO  CHECK  EQUALITY  OF  IDENT  AND  LDENT  ARRAYS  BUT  NOT  ALL 

# 

NXT=1 

ILIN1=GETWP< NXT, .REALMODE, NPPL) 

ILIN2=GETWP ( NXT , . REALMODE , NPPL) 

CORR=GETWP(NXT, . REALMODE, NLIN*NPPL) 

CORI=GETWP ( NXT, . REALMODE, NLIN*NPPL) 

# 

IF ( .OK~=OSALOC (NXT) )  GOTO  9010 

# 

CALL  COMTIN ( % 9 9 9 9 ) 

* 

CALL  FRQCTA (FDI1,FDI2, FDOl , WORK (ILIN1) ,WORK (ILIN2) , 

WORK (CORR) , WORK (CORI) , NLIN, NPPL, IEV,%9999) 

* 

CALL  PPOP 


1 


# 

RETURN 

#  ERROR  CONDITIONS 

9000  CONTINUE 

I  ILLEGAL  DATAMODE/ INCOMPATIBLE  INPUT  ARRA 

IEV=-201 2 
GOTO  9999 

# 

9010  CONTINUE 
IEV=OSGIEV ( IEV) 

GOTO  9999 

# 

9999  CONTINUE 
CALL  CLOSE (FDI1) 

CALL  CLOSE (FDI2) 

CALL  CLOSE (FDOl) 

# 

RETURN  1 
END 


•  t 


i 


# — FRQCTA  I/O  &  ALGORITHM 

# 

# IDENTIFICATION  BEEX ,  RATFOR,  VAX-II 

# 

♦PURPOSE 

#  ALGORITHM  AND  INPUT/OUTPUT  SECTION  FOR  COMMAND  THAT 

#  IMPOSES  A  FREQUENCY  DOMAIN  MEASUREMENT  CONSTRAINT 

# 

♦ENTRY  POINT 

♦  FRQCTA  ( FDI 1,  FDI2,  FDOl,  ILIN1  ,  ILIN2  ,  CORF. ,  COR  I ,  NLIN, NPPL,  IEV, 

# 

♦ARGUMENT  LISTING 

♦  FDI1  MEASUREMENTS 

♦  FDI2  OLD  RECONSTRUCTION 

♦  FDOl  NEW  RECONSTRUCTION 

♦ 

♦INCLUDE  FILES/COMMONS 

#  MACA1 

# 

♦ROUTINES  CALLED 

#  PPUSH 

#  PPOP 

#  CLOSE 

#  CPYIDR 

#  COPYDS 

#  RREAD 

#  RWRITE 

#  DSCNAM 

# 

I********************************************************* 

# 

INCLUDE  MACAl 

SUBROUTINE  FRQCTA ( FDI 1 , FDI 2 , FDOl , ILIN1 , ILIN2 , CORK, COR I , 

NLIN, NPPL, IEV, *) 

IMPLICIT  INTEGER  (A  -  Z) 

CHARACTER  FDI  1  (  .  FDLENGTH )  ,  FDI2  (  .  FDLENGTH )  ,  FDOl  (  .  FDLENG'J'H  ) 
INTEGER  IDENT(  .  IDLENGTH)  ,  LDENT  (  .  I DLENGTH )  ,  J  DENT  (  .  IDLENGTII ) 
REAL  ILIN1 (NPPL) ,  ILIN2(NPPL),  SUM,  THRESH,  SOFCOR,  DIF 
REAL  CORR( NLIN, NPPL) ,  CORI (NLIN, NPPL) 

REAL  TL,  TH,  TD 

# 

INCLUDE  TTCOM 

# 

CALL  PPUSH ( 'FRQCTA' ) 

# 

#  OPEN  INPUTFILF.S 
CALL  CPYIDR (FDI1 , IDENT, .OPNTMP, IEV, %  9000) 

CALL  CPYIDR(FDI2, LDENT, .OPNTMP, IEV,%9000) 

#  WRITE  DESCRIPTOR  RECORD  TO  TEMP  FILE 
CALL  DSCNAM( 'FRQCTA' , IEV, %9000) 

#  CALL  DSCNAM( 'FRQCTA' , IEV, %9000) 

#  SET  UP  OUTPUTFILE  PARAMETERS 

DO  1=6,19 
$( 

IF  (I  .  EQ.  6  .OR.  I  .  EQ.  7  .OR.  I  .  EQ.  13  .OR.  I  .  F.Q.  14  .OR. 
I  .EQ.  17  .OR.  I  .EQ.  19)  $( 

X=IDENT ( I ) 

JDENT(T) =X 

$) 

$) 

#  OPEN  OUTPUTFILE 
CALL  COPYDS ( FDOl , JDENT , IEV, % 9000) 


4*  4*  4* 


# 

BLKN=NLIN/2+l 

#  READ  THE  INFO  PASSKu  ALONG  IN  THE  MSMT  ARRAY 

#  LINES  1  THRU  YS  AND  HL1N-YS+2  THRU  NLIN  FOR 

#  COLS  1  THRU  XS  AND  NPPL-XS+2  THRU  NPPL 

#  CONSTITUTE  THE  MSMT  DOMAIN  IN  FREQUENCY  SPACE 
CALL  RREAD  (FDI1  ,  ILIN1 , 1 , BLKN,  IUENT,  .  WAIT,  IEV,‘i9000) 

YS=INT (ILIN1 (1) ) 

XS=INT ( ILIN1 ( 2 ) ) 

DETERMINE  MEAN  SQUARE  DEVIATION  OVER  MSMT  DOMAIN 
(EVENTUALLY  ADD  MAXIMUM  DEVIATION  OPTION  ETC) 

DO  1=1, NLIN 
$( 

DO  J=1 ,NPPL 
S  ( 

CORR(I,J) =0. 

CORI ( I , J) =0. 

$) 

$) 

SUM=0 . 

DO  BLKN=1 , YS 
$  ( 

DO  BND=1 , 2 
?  ( 

CALL  RREAD  ( FD 1 1 ,  IL INI ,  BND ,  BLKN ,  IDENT ,  .  WAIT , 1 EV , %  90  0  0 ) 

CALL  RREAD  (FDI2  ,  ILIN2  ,  BND,  BLKN,  I, DENT,  .WAIT,  IEV,%9000) 

DO  1=1, XS 
$( 

DIF=ILIN1 ( I) -ILIN2 ( I ) 

SUM=SUM+DIF**2 

IF ( BND==1 )  CORR ( BLKN , I ) =D1 F 

IF ( BND==2 )  CORI (BLKN, I) =DIF 

$) 

IF (XS>=1 ) 

$  ( 

NS=NPPL-XS+2 
DO  I=NS , NPPL 
$  ( 

DIF=ILIN1 ( I) -ILIN2 (I) 

SUM=SUM+DIF**2 

IF ( BND==1 )  CORR (BLKN, I) =DIF 

IF (BND==2)  CORI (BLKN, I) =DIF 

$) 

$) 

$) 

$) 

IF (YS>=1) 

$( 

MS=NLIN-YS+2 
DO  BLKN=MS, NLIN 
$( 

DO  BND=1 , 2 
$  ( 

CALL  RREAD(FDI1,ILIN1, BND, BLKN, IDENT, .WAIT, IEV,%9000) 

CALL  RREAD (FDI2,ILIN2, BND, BLKN, LDENT, .WAIT, IEV,%9000) 

DO  1=1, XS 
$( 


A-  I  (i 


i 


DIF  =  ILIN1 (I) -ILIN2  (I) 

SUM=SUM+DIF**2 
IF(BND==1)  CORR(BLKN, I) =DIf 
I F ( BND==2 )  COR I ( BLKN, I ) =DIF 
5) 

IF (XS>=I) 

$  ( 

DO  I=NS , NPPL 
S( 

DIF=ILIN1 ( I ) -ILIN2 (I) 

SUM=SUM+DIF**2 

IF(END=  =  1)  CORK  (BLKN,  I  )  =D1'F 

IF ( BND==2 )  CORI (BLKN, I) =DIF 

§) 

$) 

$) 

$) 

$) 

#  MEAN  SQUARE  DIFFERENCE  OVER  MSMT  DOMAIN 
SUM=SUM/ ( (2*XS-1) * (2*YS-1) ) 

#  COMPARE  WITH  TOLERABLE  MF AN  SQUARE  DIFFERENCE 

TL=0 . 0 
TH  =  1 . E+0  8 
TD=1 . 

CALL  RNGETR( 'THRESHOLD  MSV?.  1  ,  TL ,  TH  ,  TD ,  THRESH  ,  LEV ,  V.  9000) 

IF ( SUM<=THRESH)  GOTO  5000 

#  SOFT  CORRECTION 
SOFCOR=SQRT (THRESH/ SUM) 

SOFCOR® 1 . -SOFCOR 

WRITE (TTYOT, 1005)  SOFCOR 
1005  FORMAT ('  SOFCOR=  ',E9.4) 

I 

WRITE (TTYOT, 1010)  THRESH,  SUM 

1010  FORMAT ('  VARIANCE  THRESHOLD^  ’ ,E9.4,'  MS  MIFF-  '  ,E9.4) 
DO  BLKN=1,YS 
$  ( 

CALL  RREAD (FDI2,ILIN1,1, BLKN, LDENT , .WAIT, IEV,%9000) 

CALL  RREAD(FDI2, ILIN2, 2, BLKN, LDENT, .WAIT, 1EV,%9000) 

DO  1=1, XS 
$  ( 

ILINl (I) =ILIN1 (I) +SOFCOR*CORR(BLKN, I) 

ILIN2 ( I) =ILIN2 (I) +SOFCOR*CORI (BLKN, I ) 

$) 

IF (XS>=1 ) 

$  ( 

DO  I  =  NS,  NPPL, 

$  ( 

ILINl (I) =IL INI (I) +SOFCOR*CORR ( BLKN ,  I ) 

ILIN2 (I)=ILIN2 ( I) +SOFCOR*CORI (BLKN,  I) 

$) 

$) 

CALL  RWRITE (FDOl, ILINl, 1, BLKN, JDENT, .WAIT, IEV,%9000) 

CALL  RWRITE (FDOl, ILIN2, 2, BLKN, JDENT, . WAIT , IEV, % 9000 ) 

$) 

IF ( YS>=1 ) 

$( 

DO  BLKN=MS , NLIN 
$( 


A 


•  t 


CALL  RREAD (FDI2, 1LIN1 ,  1 ,  BLKN ,  LLH  NT ,  .  WAIT,  ]  IJV,¥.9000) 
CALL  RREAI)  (  FDT 2  ,  ILIN2 , 2  ,  BLKN,  I.DFNT,  .  WAIT,  IEV,%9000) 

DO  1=1, XS 
$  ( 

ILIN1  <  I )  —  TLIN1  (  I)  +SOFCOR*CORE  ( Fl.KN ,  I  } 

ILIN2 ( I) =ILIN2 (I)  +SOFCOR*CORI  (F.LKN,  I  ) 

$) 

IF  <XS>  =  1) 

?  ( 

DO  I  =  NS ,  NPPL 
$  ( 

ILIN1  (  I)  =ILIN1  (I)  +SOFCOR*COKR  ( F.LKN,  I ) 

IL IN2  ( I )  =IL  IN2  ( I )  +SOPCOR*CORI  (F.LKN,  I) 

S) 

$) 

CALL  RWRITE  ( FDOl ,  ILIN1 ,  1  ,  BLKN, JDENT,  .WAIT,  fl.V,  Y9000) 
CALL  RWRITE (FDOl , ILIN2 , 2 , BLKN, JDENT , .WAIT, IFV, V 9000) 

?) 

$) 

#  COPY  OUTSIDE  MEASUREMENT  DOMAIN 
MSM1 -MS- 1 

YSP1 = YS+1 
IF (KSM1 ' =Y  SP1 ) 

S  ( 

DO  BLKN- YS PI , MSM1 
S  ( 

CALL  RRFAD  (FDI2,  ILIN1  ,  1  ,  BLKN,  I.DENT ,  .WAIT,  I  i'M ,  V,  900  0  ) 
CALL  RRFAD  (  FD I  2 , 1 L I N2 , 2  ,  BLKN ,  I.DENT ,  .  WA 1  T ,  I  FV ,  f,  9 0 0 0 ) 
CALL  RWRITE (FDOl , IL INI , 1 , BLKN, JDFNT , .WAIT, ILV,  t 9000) 
CALL  RW R ITE  ( FDOl ,  I L I  N 2 , 2  ,  BI. K N ,  J  DENT ,  .  WA  IT,  I  F V ,  %  9 (1 0 0  ) 

$) 

S) 

GOTO  6000 

#  NO  CORRECTION,  STRAIGHT  COPY 
5000  CONTINUE 

WRITE (TTYOT, 1020)  THRESH 

1020  FORMAT  (  '  THEESIIOLD=  ',Ei!.2,'  SATISFIED') 

DO  BLKN=1 , NL IN 
$  ( 

CALL  RREAD  (FDI2  ,  ILIN1 , 1 ,  BLKN,  I.DENT,  .WAIT,  I  EV ,  Y  9  0  0  0 ) 

CALL  RREAD(FDI2,  ILIN2, 2, BLKN, I.DENT,  .WAIT,  I  EV ,  %  90  0  0  ) 

CALL  RWRITE (FDOl, ILIN1 ,1, BLKN, JDENT, .WAIT, IEV,%9000) 

CALL  RWRITE (FDOl , ILIN2 , 2 , BLKN, JDENT, . WAIT, IEV,%9000) 

# 

6000  CONTINUE 

# 

CALL  PPOP 
CALL  CLOSE (FDI1) 

CALL  CLOSE ( FDI 2 ) 

CALL  CLOSE (FDOl) 

RETURN 

#  ERROR  CONDITIONS 

9000  CONTINUE 
CALL  CLOSE ( FDI 1 ) 

CALL  CLOSE (FDI2) 

CALL  CLOSE (FDOl) 

* 


•  { 


9999  CONTINUE 
9090  RETURN  1 
END 


} 


# — MSMTSD  DRIVER 

IDENTIFICATION  BEEX,  RATFOR,  VAX-il 
PURPOSE 

DRIVER  FOR  COMMAND  THAT  GENERATES  FREQUENCY  LIMITED  MSMTS 
OVER  THE  MSMT  DOMAIN 

ENTRY  POINT  MSMTSD (WORK , ERRET) 

*****************************llt***************************x****** 

INCLUDE  MACA1 
SUBROUTINE  MSMTSD (WORK ,* ) 

IMPLICIT  INTEGER  (A  -  Z) 

INCLUDE  GIPCOM 
INCLUDE  ERROR 
INTEGER  IDENT( . IDLENGTH) 

DIMENSION  WORK (.ARB) 

# 

EQUIVALENCE  (NPPL, lDENT ( . IDNPPL) ) 

EQUIVALENCE  (NLIN, IDENT ( . IDNLINS) ) 

EQUIVALENCE  (MODE, IDENT( . IDMODE) ) 

# 

CALL  PPUSH ( ' MSMTSD ' ) 

#  OPEN  INPUTFILF 
CALL  RDKINL (FDI1 , IDENT, .OLD, IEV,%9999) 

CALL  CLOSE ( FDI 1 ) 

#  CHECK  INPUTFILF 
IF ( MODE' =. REALMODE)  GOTO  9000 

# 

NXT=1 

ILIN1=GETWP(NXT, .REALMODE, NPPL) 

ILIN2=GETWP (NXT, . REALMODE , NPPL) 

# 

IF ( .OK'=OSALOC (NXT) )  GOTO  9010 

# 

CALL  COMTIN (% 9999) 

# 

CALL  MSMTS A (FDI 1 , FDOl , FD02 , WORK ( I L INI ) ,WOKK(ILlN2) , 

NLIN, NPPL, IEV,%9999) 

# 

CALL  PPOP 

# 

RETURN 

#  ERROR  CONDITIONS 
9000  CONTINUE 

#  ILLEGAL  DATAMODE 
IEV=-2012 
GOTO  9999 

# 

9010  CONTINUE 
IEV=OSGIEV ( IEV) 

GOTO  9999 

* 

9999  CONTINUE 
CALL  CLOSE ( FDI 1) 

CALL  CLOSE (FDOl) 

CALL  CLOSE (FD02) 

# 

RETURN  1 
END 


> 


A- 20 


I 


1 


#  MSMTSA  ALGORITHM 

I  * 

'  # IDENTIFICATION  BEEX ,  RATFOR,  VAX-II 

# 

♦PURPOSE 

#  ALGORITHM  SECTION  FOR  COMMAMI)  THAT  GENERATES 

#  FREQUENCY  LIMITED  MS MTS  OVER  THE  MSMT  DOMAIN 

!  # 

;  SENTRY  POINT  MSMTSA ( FDI , FDOl , FD02 , IL INI , I L IN2 , NL IN , MPPL , I EV, * ) 

t  # 

^  ^★★ik*********************************  ******************** 

# 

INCLUDE  MACAl 

SUBROUTINE  MSMTSA ( FDI 1 , FDOl , FD02 , T L INI , I L IN2 , NL IN, NPPL , IEV, *) 
IMPLICIT  INTEGER  (A  -  Z) 

'  CHARACTER  FDI1 ( .FDLENGTH) ,  FDOl ( . FDLENGTH ) ,  FD02 ( . FDLENGTH ) 

f  INTEGER  IDENT (. IDLENGTH) ,  JDENT (. IDLENGTH ) ,  LDENT ( . IDLENGT  I! ) 

,  REAL  ILIN1 (NPPL) ,  ILIN2(NPPL),  MSV,  CONST 

f  INTEGERM  DUMMY 

\  # 

'  INCLUDE  TTCOM 


CALL  PPUSH( 'MSMTSA' ) 

#  OPEN  INPUTF1LE 
CALL  CPYIDR(FDI1 , IDENT, . OPNTMP, IEV,%9000) 
CALL  DSCNAM ( 'MSMTSA'  , IEV, % 9000) 

# 

DO  1=6,19 
$( 

IF (I  . EQ.  6  .OR.  I  .EQ.  7  .OR.  I  . EQ.  13 
I  .EQ.  17  .OR.  I  .EQ.  19)  $< 

X= IDENT ( I ) 
JDENT (I) =X 
LDENT ( I ) =X 

$) 

$) 

#  OPEN  OUTPUT!’ ILL 

CALL  COPYDS ( FDOl , JDENT, IEV,%9000) 


OR.  I  . EQ.  14 


,  OR. 


YSL  =  1 

YSH=NLIN/2-l 
YSD= (YSL+YSH) /2 

CALL  RNGETI ( 'VERTICAL  SIZE  MSMT 
DOMAIN?.  '  ,YSL,YSH,YSD,YS,  IEV,%9000)  XSI,=  .l 
XSH=NPPL/ 2-1 
XSD= (XSL+XSH) / 2 

CALL  RNGETI ( 'HORIZONTAL  SIZE  MSMT 
DOMAIN?. ' , XSL , XSH , XSD , XS , IEV, %9000)  NOIVL=0 
NOIVH=l 00000 
NOIVD=0 

CALL  RNGETI ( 'NOISE  VARIANCE?. ' , NOIVL, NOIVH , NOIVD, NOIV, IEV,%9000) 

# 

CONST=SQRT ( 12 . *NOIV) 

IF ( NOIV==NOIVD)  CONST=0 . 0 
DUMMY=345  89 
DO  BLKN=1 , NLIN 
$( 

CALL  RREAD (FDIl , ILIN1 , 1 , BLKN, IDENT, .WAIT, TEV,%9000) 

DO  1=1, NPPL 
$( 

ILIN1 (I) =ILIN1 ( I) +CONST* (RAN (DUMMY) -. 5) 

$)  A- 2  I 


3 


CALL  RWRITE (FD01 , ILINl , 1 , BLKN, JDKNT, .WAIT, ] EV,%9000) 
CALL  RREAD (FDI1 , ILINl ,2 , BLKN, IDENT, . WAIT , IEV, %  9000 ) 
DO  1=1, NPPL 
$  ( 

ILINl  (I)  =ILIN1  (I)  +CONST*  ( RAN  (DUMMY  >-.  r> ) 

5) 

CALL  RWRITE  (FD01. , ILINl , 2 , BLKN, J DENT,  .WAIT,  IEV,%9000) 

$) 

# 

BLKN=NLIN/2+l 
DO  1=1 , NPPL 
$< 

ILINl (I) =0. 

$) 

ILINl(l)=YS+.02 
ILINl (2)=XS+. 02 

CALL  RWRITE (FDOl, ILINl, 1, BLKN, JDENT, .WAIT, 1EV,%9000) 

# 

CALL  CLOSE (FDOl) 

CALL  CLOSE (FDI1) 

CALL  CPYIDR(FDI1 , IDENT, . OPNTMP, IEV,%9000) 

CALL  COPYDS (FD02 , LDENT, IEV,%9000) 

XSP1=XS+1 

NS=NPPL-XS+2 

NSM1 =NS~1 

YSPl=YS+l 

MS=NLIN-YS+2 

MSM1=MS-1 

# 

MSV=0 . 0 
DO  BLKN=1 , YS 
$( 

CALL  RREAD (FDI1, ILINl , 1 , BLKN, IDENT, .WAIT, 1EV,%9000) 
CALL  RREAD  (  FDI 1 ,  II,  IN2 , 2  ,  BLKN ,  IDENT,  .WAIT,  IEV,  V.  9000 ) 
DO  1=1, XS 
$  ( 

MS  V=MSV+ ILINl  ( I)  *  ILINl  (I)  +ILIN2  ( I)  *II,TN2  (  I) 

S) 

DO  I=NS , NPPL 
$( 

MSV=MSV+ ILINl (I) * ILINl ( I) +ILIN2 ( 1) *1L LN2 ( I ) 

$) 

DO  I=XSP1 , NSM1 
$  ( 

ILINl ( I ) =0 . 0 
ILIN2 ( I ) =0 . 0 

$) 

CALL  RWRITE (FD02, ILINl, 1, BLKN, LDENT, .WAIT, IEV,%9000) 
CALL  RWRITE (FD02, ILIN2 ,2, BLKN, LDENT, .WAIT, IEV, % 9000) 

$) 

DO  BLKN=MS , NLIN 
§( 

CALL  RREAD(FDI1, ILINl, 1, BLKN, IDENT, .WAIT, IEV,%9000) 
CALL  RREAD (FDI 1 , ILIN2 , 2 , BLKN, IDENT, .WAIT, IEV, %  9000) 
DO  1=1, XS 
$( 


S) 


a-  t: 

_ A 


MSV=MSV+ILINl ( I ) * ILINl ( I ) +ILIN2 ( I ) * 1LIN2 ( I) 


DO  I=NS , NPPL 
$( 

MSV=MSV+ILIN1 (I) *ILIN1 ( I ) +ILIN2 ( I ) * IL 1 N2 ( I ) 

$) 

DO  I=XSP1 , NSM1 
$  ( 

ILIN1 ( I ) =0 . 0 
ILIN2 ( I ) =0 . 0 

$) 

CALL  RWRITE  (FD02  ,ILIN1,1,  BLKN,  LDENT,  .WAIT,  IEV,%9000) 
CALL  RWRITE <FD02,ILIN2, 2, BLKN, LDENT, .WAIT, IEV,%9000) 

$) 

DO  BLKN=YSP1 , MSM1 
$( 

DO  1=1, NPPL 
$( 

ILIN1 ( I ) =0 . 0 
ILIN2 ( I ) =0 . 0 

$) 

CALL  RWRITE (FD02, ILIN1 , 1 , BLKN, LDENT, .WAIT, 1EV,%9000) 
CALL  RWRITE (FD02,ILIN2, 2, BLKN, LDENT, .WAIT, IEV,%9000) 

$) 

MSV=MSV/( ( 2*XS-1 ) * (2*YS-1)  ) 

WRITE (TTYOT ,1010)  MSV 

1010  FORMAT ( 1  MEAN  SQUARE  VALUE  OF  MSMTS=  ’ ,E10.4) 

# 

CALL  PPOP 
CALL  CLOSE (FDI1) 

CALL  CLOSE (FDOl) 

CALL  CLOSE (FD02) 

RETURN 

#  ERROR  CONDITIONS 

9000  CONTINUE 
CALL  CLOSE (FDI1) 

CALL  CLOSE (FDOl) 

CALL  CLOSE (FD02) 

# 

9999  CONTINUE 
9090  RETURN  1 
END 


A-l’ 
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